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We study a lattice model for the spreading of fluid films, which are a few molecular layers thick, in 
narrow channels with inert lateral walls. We focus on systems connected to two particle reservoirs at 
different chemical potentials, considering an attractive substrate potential at the bottom, confining 
side walls, and hard-core repulsive fluid-fluid interactions. Using kinetic Monte Carlo simulations 
we find a diffusive behavior. The corresponding diffusion coefficient depends on the density and is 
bounded from below by the free one-dimensional diffusion coefficient, valid for an inert bottom wall. 
These numerical results are rationalized within the corresponding continuum limit. 
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I. INTRODUCTION 

In recent years, substantial progress has been made 
in the development of the "lab on a chip" concept, i.e., 
the integration of many physical and chemical processes 
(e.g., transport through micro-channels, mixing of differ- 
ent fluids, chemical reactions) into a single device; en- 
tire laboratory setups, like a gas chromatograph, have 
been miniaturized on a single chip (for a review see, e.g., 
Ref. [l|). In this context, microfluidics is becoming a 
standard tool in many applications, ranging from biology 
(see, e.g., Ref. 0) to the handling of toxic or rare sub- 
stances. Further scaling down to nanofluidics is expected 
to take place in the future [3j. Already now it is possi- 
ble to sculpture channels with lateral dimensions of few 
tens of nanometers Ij] (for a review on such fabrication 
processes see Ref. @) and carbon nanotubes have been 
proposed as possible pipes in nanofluidics [1, 0] • Chemi- 
cally patterned substrates have also been suggested as a 
solution for directed transport, gating, mixing, or sepa- 
ration of fluids at the micro- and nano-scale [8j. In this 
case the channel consists of a strip of wettable material 
embedded in a non-wettable substrate so that the fluid 
flows along the wettable region and is laterally confined 
by the chemical contrast. 

If one of the dimensions of a fluid film is comparable to 
the size of the fluid molecules, a hydrodynamical descrip- 
tion of the film is no longer justified 0, [l(| El- In this 
case the discrete nature of the fluid becomes important 
and the fluid cannot be treated as a continuum in the con- 
fined direction. In order to investigate such systems one 
possible approach is to carry out computer simulation of 
discrete models, e.g., molecular dynamics, kinetic Monte 
Carlo (KMC), or lattice Boltzmann simulations; recent 
work in this direction includes fluids in carbon nanotubes 
[(| and on chemically patterned substrates [l2l ]. 

With the scaling down of microfluidic devices one has 
to deal with and may exploit the ultrathin precursor 
film which spreads ahead of the bulk fluid. Experi- 
mental studies have shown that in some cases such pre- 
cursor films have molecular thickness [H, 0, EH, Gil 



UtI [H, 0, The spreading of such monolayers has 
been studied using a two-dimensional lattice gas Ising 
model [ToL |2ll |22| in which a half-space is occupied by a 
particle reservoir. Recent KMC simulations and a con- 
tinuum analysis [23[ of that model provided results in 
good qualitative agreement with available experimental 
data, and a further extension to the case of chemically 
patterned substrate has been proposed (24|. 

Fluids in narrow channels have been investigated the- 
oretically in the context of single-file diffusion, i.e., when 
fluid particles cannot overtake each other (see, e.g., 
Refs. |i, |H H3, [H, [H, [H SH HI H 11 ) . Such systems 
show the interesting feature of non-diffusive behavior of 
tracer particles, which stimulated experimental (see, e.g., 
Refs. (Si, [H) and numerical [H IM [H HH interest. 
Here we present a lattice model for ultrathin films in 
which multiple occupancy of a site is allowed ( gene raliz- 
ing the single-occupancy model of Refs. [ToL [2lTT23| ^ and 
in which the substrate-particle attractive interaction is 
decaying as a power law, whereas the particle-particle in- 
teraction is assumed to be hard-core repulsive only. This 
mimics the case in which the fluid-substrate interaction 
strongly dominates over the actual attractive long-range 
part of the fluid-fluid interaction. Based on the phenom- 
ena occurring in this minimalistic model, the extension 
to the case in which the attractive part of the fluid-fluid 
interaction is relevant will be presented elsewhere. We 
shall restrict our analysis to a one-dimensional model, 
which can effectively describe fluids in extremely narrow 
channels with a width which is less than twice the parti- 
cle diameter. The sidewalk act to confine the particles. 
The corrugation of the substrate potential both at the 
bottom and at the sides is incorporated effectively by 
considering a lattice model for the particles. Due to the 
small thickness of the channel the transversal variation 
of the substrate potential can be ignored. This model 
is supposed to mimic not only molecular fluids but also 
colloidal particles in solution, with the colloidal particle 
setting the length scale. 

We discuss both the initial dynamics, in which a fluid 
film fed by a reservoir gradually fills the channel, and the 
steady state, in which the fluid film in the channel is in 
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contact with two reservoirs at different chemical poten- 
tials. 

The paper is organized such that in Sec. [TT] we define 
the model whereas in Sec. IIIII the results of our Monte 
Carlo simulations are presented. The analyses of the 
diffusion-like dynamics and of the steady-state proper- 
ties are presented in Sec. IIV1 In Sec. [V] we discuss the 
mean-field continuum limit of the model and rationalize, 
within this approximation, the results for the diffusion 
coefficient presented in Sec. |IV] Section |Vl] summarizes 
the main findings and provides our conclusions. 

II. THE MODEL 

Before specifying the rules defining the model, we fur- 
ther describe the general physical picture of the type of 
systems we have in mind. 

As stated above, the fluid is assumed to be confined 
to a narrow, effectively one-dimensional channel. The 
sidewalls are very high compared with the fluid particle 
diameter, so that the fluid cannot spill out of the chan- 
nel. The channel walls act on the particles such that only 
the vertical variation of the substrate potential matters. 
The left and the right end of the channel are connected 
to a feeding and absorbing particle reservoir, respectively, 
and the channel is initially empty. The fluid film inside 
the channel is taken to be compact, i.e., molecules are 
densely packed to form vertical columns without vacan- 
cies. This corresponds to the case in which the substrate 
is strongly attractive and vacancies inside columns are 
eliminated on a time scale much shorter than the typical 
time for exchanges of particles between columns. We 
describe these exchanges in terms of rates, which are 
related to the change of the energy of the system due 
to the corresponding move. Particle exchanges between 
columns and particle insertions and removals near the 
reservoirs are the only processes we consider. We assume 
that neither evaporation nor condensation takes place in- 
side the channel. This minimalistic model aims at identi- 
fying general aspects and the main qualitative features of 
fluids spreading in a strongly confined geometry, rather 
than providing an accurate description of a particular 
physical situation. 

Inspired by this picture, in Subsec. Ill Al we specify the 
configuration space and the corresponding energy func- 
tion. In Subsec. Ill Bl the dynamical rules (i.e., the al- 
lowed changes of the configurations and the associated 
rates) governing the time evolution in the bulk are dis- 
cussed, whereas Subsec. Ill CI deals with the definition of 
the dynamics at the feeding and absorbing boundaries. 

A. Configurations and Hamiltonian 

The model is defined on a one-dimensional (D = 1) 
lattice, with sites [0, . . . , (L + l)a]. The distance a be- 
tween two consecutive sites is assumed to be equal to the 
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FIG. 1: A typical configuration of the model. The possi- 
ble moves in the bulk are indicated by straight arrows, while 
the curved arrows denote reservoirs-system exchanges. The 
substrate, including the exclusion zone of its top layer, cor- 
responds to the hatched area. The grey areas at x < and 
x > (L + l)a indicate reservoirs and the fluid particles are 
shown as circles. 

effective diameter of a fluid particle, which is set by the 
hard-core repulsion between the particles. In the follow- 
ing, the site indices and the distances will be expressed 
in units of a. The two sites indexed with and L + l are 
the boundaries of the left and right particle reservoirs, re- 
spectively. The other sites, [1, . . . ,L] (called 'bulk' in the 
following) , represent the channel of length L (assumed to 
be long, i.e., L^> 1). 

At every site x € [0, L + I] an occupation number 
n x G No specifies the number of particles piled in the 
column x (see Fig. Q]). Within the column, particles cen- 
ters are located at integer y positions. Accounting for 
fluid-substrate hard-core repulsion we consider the posi- 
tion y = as passing through the centers of the particles 
forming the top layer of the substrate. If the diameter 
of the substrate particles differs from that of the fluid 
particles, one may introduce an extra parameter to char- 
acterize the position of the fluid-substrate contact layer; 
for simplicity, however, we assume that the fluid particles 
in the first layer are located at y = 1 (see Fig. [T]). 

The substrate is assumed to be uniform, and, consis- 
tent with our one-dimensional model, two-dimensional 
semi-infinite in the y < - direction. We denote the 
(attractive) pair interaction between a substrate particle 
(p) 

and a fluid particle by U s j , resembling dispersion forces: 

where d is the dimensionless distance in units of a be- 
tween the substrate particle located at (a/, —y' < 0), and 
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the fluid particle located at (1 < x < L, y > 1). In the 
case of pairwise additive interactions, for a semi-infinite 
substrate (x' £ R, y' £ R + ) and in the continuum limit 
(d 3> 1), this leads to a total substrate potential 



U 8 f{y) = -w s f / dy / dx 
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where w' s j = ^jw s f. Within this ansatz, the particle- 
substrate interaction in Eq. ([3]) depends on the height y of 
the particle only. The energy of the fluid configuration 
{ni,...,ni,} exposed to the substrate potential U s f is 
thus given by 



(4) 



where the inner sum is defined to be if n x = 0. 
Note that, following the discussion at the beginning of 
the present section, we assume that columns are always 
densely packed, so that configurations are as depicted in 
Fig. [TJ since the configurations are characterized com- 
pletely by a succession of numbers {n\, . . . , ul}, the en- 
ergy is a function of these numbers only, as in Eq. ^ . 

The same form [Eq. (pj] of the pair potential is as- 
sumed for the fluid particle - fluid particle interaction, 
where the corresponding interaction strength is denoted 
by Wff. Each pair of particles separated by a distance 
dff > 1 contributes to the particle-particle energy, so 
that the total energy due to particle-particle interactions 
can be written as 



^ L L n x n x i 

ff = 2 ^ ^ S £ 

i=li'=lV»=l» 1 '=l 

# (dff = ^{x-x>y + {y x -y*>) 2 ) 



(5) 



B. The rates and the dynamics 

In this subsection we define and discuss the rates which 
govern the stochastic dynamics in the bulk, i.e., for x g 
[1,L]. The dynamics at the boundaries, x = and x = 
L + 1, will be discussed in the following subsection. 

We assume that each particle in the column x may 
jump into one of the nearest neighbor (NN) columns x+l 
or x — 1. We introduce the rate J'cc" (?A 2/')j which is the 
rate for a particle in column x and at given height y, to 
jump to the next column x + l and at height y' . Within 
our aforementioned model assumption this process in- 
volves an instantaneous column height reduction by one 
in column x and a column height increase in column x+l. 
This also means that the jumping particle is considered 
to be able to squeeze into column x + l at position y' by 
pushing the particles above this position up by one unit; 
on the other hand if y' is above the top particle of col- 
umn at + 1, it falls down in order to form again a compact 
column. Accordingly the configurations C,C, 



C = {ni, . . . ,n x ,n x+ i, . . . ,n L } 
C = {ni, ...,n x - l,n x+ i + 1, . . . ,n L }, 
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represent the initial and the final configurations, for any 
pair y,y' . Analogous considerations can be carried out for 
moves from x + 1 to x, where the above configurations 
are interchanged. Therefore, within our model, the rates 
fCC'{y-,y ! ) depend only on the initial and final configu- 
rations C and C", respectively. Accordingly, the depen- 
dence on y,y' is dropped. We introduce the dimensionless 
rate ucc, which is also assumed to depend only on C,C, 



UCC 
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where v§ fixes the time-scale of the model and we assume 
it to be independent of the source and target columns fill- 
ing, i.e., the same for any particle in the source column. 
(vq can be interpreted as the rate for a particle to jump to 
the NN column if the energy happens to be unchanged 
by the move.) In the following, times are measured in 
units of vq, i.e., the dimensionless simulation time t cor- 
responds to an actual time t a = t/vQ. 

We choose the rates u for all possible moves from col- 
umn x to column x + l such that detailed balance 



with uf]{Q) = and the sums over y x and y x i are taken 
to be zero if n x = or n x ' — 0. The total energy function 
is H[C] = Hff + H s f where C = {ni, . . . , n^} character- 
izes completely each configuration. Note that this part 
of the Hamiltonian is restricted to the bulk; in general 
the reservoir-bulk interactions should also be accounted 
for separately. In the special case of the absence of long- 
range particle-particle interaction, i.e., Wff =0, both the 
bulk and the reservoir-bulk contributions vanish, and the 
energy is H [C] — H s j. As mentioned in the Introduction, 
in the following we shall discuss only this situation; the 
case Wff ^ will be presented elsewhere. 
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is obeyed, where AH[C,C] = H[C] - H[C] is the en- 
ergy difference between the final (C) and the initial (C) 
configuration. Detailed balance has been chosen in order 
to ensure that thermal equilibrium is reached in the long- 
time limit, if the two reservoirs of particles at the right 
and the left of the channel are set to the same chemical 
potential. A possible choice that satisfies the detailed 
balance condition is 
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The chosen form of the rates [Eq. |9|)] includes both 
"slow" (AH > 0) and "fast" (AH < 0) processes, and 
we implicitly assume that it captures essential features 
of the real dynamics. 

The rate iicc is the same for any particle in the source 
column x, so that the total rate ucc for a column to 
decrease its occupation number by one, while a given 
NN column increases its own occupation number by one, 
is 



n x u C c 
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The rates in Eq. (fTQ|) are defined on the space of config- 
urations specified by occupation numbers only. Detailed 
balance still holds and the corresponding Boltzmann sta- 
tistical weight is 



p B (ni, ...,n x ,...,n L )oc 



-0H 
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which accounts for "particle undistinctness" by dividing 
the Boltzmann factor by n\ !,..., nx! where n x \ is the 
number of choices to label the n x particles in each x € 
[1, . . . , L] column. In the case of the particle-substrate 
interaction described by the Hamiltonian in Eq. the 
rate in Eq. (flQ|) has the following explicit form: 



u (n x 



_i) = n x exp 
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This formula emphasizes that u depends only on the oc- 
cupation numbers of the initial and the target column. 
The notation is such that the first argument stands for 
the source column (here located at x with occupation n x ), 
while the second argument represents the target column 
(here at x + 1 with occupation n x+ i). 

Assuming that the dynamics leads to a diffusion-like 
behavior (as will be discussed in Sec. IIV[) , some quali- 
tative features of the diffusion coefficient as a function 
of the local density can be anticipated from the general 
properties of the rates in Eq. (fT2| . First, consider the 
situation in which both n x = n and n x ±\ — m are large 
compared to (fiw 1 ^) 1 ^ . Then the exponent in Eq. (fT2|) is 
very small and u(n,m) — ► n, so that the model reduces 
to free particles diffusing in D = 1. The same conclu- 
sion holds for n = m + 1, in which case the exponent 
is zero leading to u(n,m) = n. In general, u(n,m) ^ n 
if n ^ m + 1; accordingly, jumps from high columns to 
low columns are faster than in the free case, while the 
opposite processes are slower. This means that diffusion, 
which tends to smooth out density gradients, is enhanced 
by the exponential factor in Eq. (|12|). Since at low den- 
sities most of the configurations are composed either of 
empty columns or of columns occupied by one particle, 
the most probable rate is w(l,0) — 1, which results in 
free diffusion at low densities. 

These considerations lead to the conclusion that the 
diffusion coefficient is expected to exhibit a peak at rel- 
atively low densities, because the rates exhibit the max- 
imal difference with free diffusion rates if the target col- 
umn is empty. 



Before passing to the definition of the dynamics at the 
boundaries, we briefly comment on similar models which 
have been considered in the literature. In Refs. [35|,[36| a 
class of dynamical models, to which our model belongs, is 
introduced and studied. In this class of models the rates 
depend on both the source and the target column, they 
do not necessarily satisfy detailed balance, and jumps 
occur not only between NN. (These models are known 
in the literature as misanthropic processes.) The main 
result of Refs. [35L [36j is that under certain conditions in 
the infinite square lattice it is possible to obtain an exact 
expression for the steady-state distribution. A concise 
summary of these results can be found in Refs. [13, HU, 
where their relevance for the non-equilibrium dynamics 
of interacting particles has been stressed. Applied to our 
case, the results in Refs. [35L l36j recover the equilibrium 
Boltzmann distribution in Eq. (jlip , with the Hamiltonian 
defined in Eq. but do not provide information on 
the dynamics and steady-state distribution if chemical 
drive, caused by different chemical potential for the two 
reservoirs at the boundaries, is applied. 



C. Dynamics at the boundaries 

We consider now the dynamics at the boundaries and 
discuss two possible implementations. The first choice 
is to fix the occupation number of the columns and 
L + 1 at values no and nx+i, respectively, and to impose, 
with some additional assumptions, the same dynamics 
as in the "bulk" . The boundary dynamics changes the 
occupation number ri\ of the first column of the system, 
according to Eq. (|12p , while the occupation number no 
is unchanged by the move, and the same holds at x = L. 
One can physically motivate such a choice by assuming 
that the particle exchanges within the reservoir are so 
fast, so that a particle extracted from the reservoir is 
immediately replaced. Under this assumption the density 
of particles in the reservoir is simply no. While this choice 
seems to be rather natural, as explained in the following 
the equilibrium (i.e., for no = n^ + i) properties display 
an unexpected feature, i.e., a jump discontinuity in the 
density between the reservoirs and the system. 

The total Hamiltonian in Eq. is a sum of single- 
column terms, so that the equilibrium grand canonical 
distribution factorizes: 

1 L 

P eq ({rn, . . . , n L }) = — TT Pw (n k ) e-» n * , (13) 

where 



Z(w,n)= {J2p w (n)e-^ n 



(14) 



Vn=0 



is the total partition function, fi — (3jl is the dimension- 
less chemical potential, with ft, as the actual chemical po- 
tential, p w is a non-normalized single-column statistical 
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weight, 



p w {m) = — - exp [2wh(m)], 



w = (3w'/2 is a dimcnsionless quantity, and 



1 



0, for n = 0. 



(15) 



(16) 



The chemical potential n controls the mean density of 
particles in the system, 
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Peq(z) = -J 1 = ~d z In [Z (w, z)} 



^2 np w (n)z r 



n>0 



p(n)z r ' 



where N eq is the mean total number of particles in the 
system, and z = e _Al is the fugacity. Detailed balance for 
the rates at the left reservoir reads 



Peq({ni 



,ni})u(no,ni) 



P eq {{ni 



1. 



«(ni + l,n ) 



(18) 



where the probability per unit time of inserting a par- 
ticle into the column at x = 1 is compared with the 
corresponding probability per unit time of removing the 
particle in the same column. A similar condition has to 
hold at the right end x = L of the system. Combining 
Eq. HU with Eqs. JUJ) and (13]) leads to 



no exp 
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no K + I) 4 



(19) 



Equations (fT7|) and (|19|) give the equilibrium density p eg 
in the system as a function of uq. As expected, if no is 
large p eq coincides with no: in Eq. (|19p no 3> w 1 / 4 im- 
plies z ~ no, and in Eq. (fT7|) no ^> w 1 ^ 3 implies p w {ri) »s 
exp (2w£(4))/n!, so that for n 3> max (u; 1 / 3 , if 1 / 4 ) 
Eq. (JTTJ) reduces to 



d 



Peq 



z — e 

az 



■ w no, no > max(w 1,,3 ,M) 1 ' 4 ). 

(20) 

In the range of substrate potential strength we investi- 
gated (0.5 < w < 5) the approximation p eg no is 
valid if n > 5. As the substrate potential strength is 
increased, this threshold increases, while it tends to zero 
for w — > 0. For densities lower than the threshold, the 
reservoir occupation number no does not coincide with 
the density p eq of the equilibrium system as obtained 
from Eqs. fB] l -([T7 |) and (0 (see Fig.©. For example, 
in the case no = 3 one has p eq ~ 3.19 for w = 2 and 
p eq ~ 3.23 for w — 4. These densities are in very good 
agreement with the simulation data for the density p\ in 
the first column (see, c.f., Fig. GUb)). In the simulations 
we investigated a non-equilibrium situation in which the 



FIG. 2: The equilibrium density p eq [Eqs. lfT3 ]l -[[T7 |) and ([19]) ]. 
divided by the reservoir occupation number no, as a function 
of no for uu = 2 and w = 4, respectively. 



two reservoirs at the boundaries have different occupa- 
tion numbers (n ^ n-i+i); nevertheless we recover the 
equilibrium density in the first column. This shows how 
it is possible to control the densities at the first {x = 1) 
and last (x = L) site by varying the occupation numbers 
no and ul+i- In order to obtain arbitrary densities, it 
is necessary to take no and nj, to be continuous, thus 
loosing the direct physical interpretation of these param- 
eters. Moreover, as shown in Fig. [2 p eq drops sharply 
for n < 0.8, and in the range < p eq < 0.8 a high 
numerical accuracy would be required to determine the 
corresponding value of no . 

These two problems can be solved by generalizing the 
dynamics at the boundaries as follows. In Eq. (TT2"]) the 
terms depending on no and Ul+i, i.e., the properties of 
the reservoirs, are replaced by constants a, 7, 6, and k 
in the following way: 



u a (n{) = a exp 
us{n L ) = S exp 



(m + 1) 4 

w 



(n L + l) 4 



, n 7 (ni) = 7 ni exp 
, u K (n L ) — nn L exp 



w 

< 
w 



(21) 



where u a is the rate for particle insertion into column m 
from the left reservoir, u 7 is the rate for particle removal 
from column n\ into the left reservoir and us,u K are the 
corresponding rates at x = L. Imposing detailed balance 
[Eq. ([15])] for the rates defined in Eq. ([IT]) gives 



7 K 



(22) 



In the non-equilibrium case the fugacity of the right reser- 
voir, denoted as Zl+i — S/k, and the one of the left 
reservoir, i.e., zq — a/7, are different (zq 7^ Zl+i). Us- 
ing these two fugacities in Eq. (fTT]) . the densities of the 
two corresponding equilibrium systems are found; we de- 
fine them to be the reservoir densities. In simulations 
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we proceed backwards: first we choose po = Peq(zo) and 
Pl+i = Peq( z L+i), and then we find the corresponding 
ratios by inverting Eq. (|17|) . Setting 7 = k = 1 implies 
a = zq, S = zl+i so that the inversion is simpler. 



III. 



MONTE-CARLO SIMULATIONS 



The continuous time dynamics defined by the rules 
described in Subsecs. Ill Bl and III CI is simulated using a 
Kinetic Monte Carlo (KMC) method [39] . At every step 
an increment At for the time variable is drawn from the 
distribution 



P(At) 



1 



S (n , . . . , n L+ i) 



exp[5 (n , • . . ,n L +i) At], 



S (n , . . . , n L+ i) = 2J [u(n x ,n x+ i) + u(n x+ i, 



(23) 



where S is the total rate to leave the configuration 
{no, . . . , ul + i}. The move to perform is then chosen ac- 
cording to the weight u/S of its rate. We used a classical 
N-fold way algorithm [4(| , which has the advantage that 
the selected moves are accepted without rejection. The 
model depends on four parameters: the substrate interac- 
tion strength w = w' sp (3/2 in units of the thermal energy, 
the two boundary densities no and n^ + i, and the length 
L of the system. The simulations have been performed 
up to a maximum time r tot and quantities have been 
measured after the initial time To. The simulations cov- 
ered both the spreading and the steady-state regime and 
in both cases we sampled the same set of quantities. In 
order to keep the notation simple we indicate averages 
always with (•), but, as described in the following, the 
meaning of the symbol is different in the two situations 
considered. 

The mean density of particles, i.e., the density at site 
x and time t is defined as 



p{x,t) = (n x (t)), 



(24) 



while the total mean number of particles (or mean total 
mass) is 



M(t) = £> x (t)>. 



(25) 



The transport properties have been studied using the in- 
tegrated particle current at site x and time t, 

J{x, t, At) = An x , x+1 (t, t + At)- An x+hx (t, t + At), 

(26) 

where An„/ (t, t+At) is the number of particles jumping 
from site x to site x' within the time interval At. We 
define a mean instantaneous current as 



(j(x,t))= A lim 



(J{x,t,At)) 
At 



(27) 



A. Spreading 

For the spreading regime the right reservoir has been 
converted into a particle sink by setting ul+i = 0, while 
no = 11 and L — 1000. In the simulations performed 
with these values of the parameters the leakage of parti- 
cles through the sink is negligible for times t < 10 4 . We 
studied both the initial-time dynamics by setting to = 
and Ttot < 10 4 and the long-time behavior in which both 
reservoirs play a role (see, c.f., Fig [4] w = 0.5); in this 
latter case we have chosen r tot — 5 x 10 4 and t = 10 4 
in order to reduce the CPU and memory requirements. 
In the spreading regime we implemented the simulation 
average by drawing different sequences of (pseudo-) ran- 
dom numbers while keeping the initial condition fixed so 
that in this case (•) is the ensemble average; the typical 
number of runs we averaged over is 2000. We studied the 
shape of the density profile p(x,t), defined in Eq. (I24|) . 
as a function of the interaction strength w in the range 
0-5 < w < 1.5. 

Since for w = the dynamics reduces to free diffusion, 
it is natural to check if in the general case (i.e., w 7^ 0) 
the profiles show a diffusive scaling. Plotting them as 
a function of A = xj\ft we indeed obtain a collapse of 
data measured at different times, as shown in Fig. [3] 
The rescaled profiles are similar to the one for free diffu- 
sion, except for a small bend at densities around 1 (see 
FigE^b)), which depends on the interaction strength w. 
The diffusive scaling is confirmed by the time evolution 
of the total mass [Eq. (f2"5)) ] shown in Fig. 0] which is ex- 
pected to evolve as M oc \fi. We observe deviations from 
this behavior only at short times, when the boundary dy- 
namics dominates, and at long times, when the leakage 
of particles through the sink becomes relevant. 



B. Steady state 

The steady state is reached after running the simula- 
tion for an initial thermalization time To (to ~ 10 5 for 
L = 1000, chosen by checking that for t > To the ob- 
scrvables are time independent) and saving the configu- 
rations generated every sampling time interval t s , with 
t s = 200 for L — 1000. The average (•) for the observ- 
ables defined above is taken over this set of configura- 
tions. The choice for t s is a compromise between speed 
and having as small correlations between the N s measure- 
ments as possible. We assume that the total simulation 
time T to t — N s t s + tq is sufficient to explore a significant 
part of the phase space, so that the performed average 
coincides with the average over the (unknown) steady 
state distribution. Note that this assumption is justified 
because no signs of dynamical phase transitions (which 
would introduce extremely long time scales) are found in 
the simulations. 

The density p(x) [Eq. (f2~4")) ] exhibits a profile smoothly 
interpolating between the two reservoirs (see Fig. [5]) and 
slightly deviating from the corresponding free diffusion 



7 




FIG. 3: Time- dependent density profiles for spreading in a system with L — 1000, for no = 11, tiL+i = 0, and w = 1.5 at times 
t = 5 x 10 3 (+), 10 4 (□), 2 x 10 4 (■), and 2.5 x 10 4 (0) as a function of x (a), and of A = x/y/t (b), respectively. The solid 
line indicates the scaling function for free diffusion (i.e., w = 0). 
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profile which is a straight line. This deviation is consid- 
ered in more detail in Appendix [B] where its dependence 
on both the interaction and the reservoirs densities is an- 
alyzed. In the steady state both reservoirs play a role and 
finite-size effects have to be checked; it turns out that for 
L > 200 there is no detectable dependence of the data 
on the particular value of L other than a trivial rescaling 
of the density profile. 

We have also determined the current (j) defined in 
Eq. (|2T|) : in the steady state (j) does not depend on t, 
so that J(x, t, At) — (j) At for any sufficiently large time 
interval At. J can be obtained by measuring the flux of 
particles between any two sites x and x + 1, because in 
the steady state the current {j) does not depend on x due 



to local particle conservation. The instantaneous steady 
state current (j{x)) (a dependence on x is indicated to 
recall the random fluctuations around the mean value 
(j)) can then be obtained from 



(i(z)> - 



J(x,0,-r tot ) - J(x,0,tq) 
N x t, 



(28) 



Note that in Eq. I|28p the current integrated over the 
thermalization time, J(x, 0, To), which depends on x and 
t, has been subtracted. We have calculated (j(x)) for 
x G [1, . . . , L] leading to (j) = L' 1 £x=i0'(*)>- 



IV. ANALYSIS OF THE SIMULATION DATA 
A. Methods to determine the diffusion coefficient 

Guided by the results of the KMC simulations of the 
microscopic model we expect that a continuum (in space 
and time) description for the behavior of the model at 
long times and large spatial scales is possible, i.e., that 
the hydrodynamic limit exists and is well defined. A 
rigorous proof has been provided for a small number of 
models (see, e.g., Ref. [4l|). In the present case such an 
explicit derivation appears to be a difficult task. 

Assuming that the particle density p and the current 
(j), defined in Eqs. (f2"4")l and (|27[) . are smooth functions 
of the position x and of the time t, the local conservation 
of particle density in the bulk (which is implicit in the 
dynamics of the model) is expected to take the form of a 
continuity equation: 



d t p(x,t) = -d x (j(x,t)). 



(29) 



The results of the simulations strongly indicate a dif- 
fusive scaling at long times and large spatial scales, i.e., 
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FIG. 5: Steady-state density profiles in a system of length L — 1000 for (a) w = 0.4 (□) and w = 1.5 (0), no = 8, ul+i = 0; 
(b) w = 2 (□) and w — 4 (©), no = 3, ul+i = 0. The inset in (b) is a close-up view of the vicinity of the left reservoir showing 
p eq ~ pi > no [p eq ~ 3.19 for w = 2 and p eq m 3.23 for w — 4; see Eqs. (fT7|) and (fl"9|) and the discussion in the main text]. In 
both (a) and (b) the free diffusion is indicated by a full line. 



p(x, t) = p(x/y/t), suggesting that the dynamics amounts 
to non-linear diffusion: 



(j(x,t)) = -D(p)d x p{x,t) => 
d t p(x,t) = d x [D(p)d x p(x,t)} 



(30) 



Based on the density profile p(x, t) from the simulations, 
it is possible to extract the function D(p) from the data in 
the spreading and the steady state regime, respectively, 
as follows. 

• Spreading. Using the scaling behavior p(x, t) — 
p(X = xfy/i) (see the results in Sec. [ITTJ), Eq. ([30]) 
reduces to 



D(p{X))—p(X) 



(31) 



Assuming that ^D(p(X)) and p(X) vanish for A — > 
oo, which is supported by the simulation data 
(L > 1 is an approximation for L — > oo), inte- 
grating Eq. (|31|) and inverting p(X) into A(p), one 
finds 



D(p) = ±X(p) [ P dpX(p). 
dp 7o 



(32) 



This method might be inaccurate for small densi- 
ties due to a systematic effect. For x « L the profile 
bends in order to fulfill the condition n^+i = 0, so 
that its derivative -^p is larger than the derivative 
of a profile poo in the infinite lattice: gjPoo < j\P- 
The integral in Eq. (|32p is also underestimated, 
since X(p — > 0) oo. These effects lead to an 



underestimated diffusion coefficient for small val- 
ues of p. The most severe effect is probably due 
to the derivative -^p(X), while the errors in the 
integral can be partially corrected by using larger 
lattice sizes. 



• Steady state. In the stationary state the current 
and the density are constant with respect to time, 
so that Eq. $5U§ leads to 



D(p) 



(33) 



where p'(x) = 4^p(x). Note that the computation 
of D in this regime does not require any assump- 
tions on p and p ', as in the previous case, and 
therefore no systematic deviations from the actual 
diffusion coefficient are expected. 



• Steady state in quasi-equilibrium. A steady 
state deviating only slightly from the equilibrium is 
realized by imposing reservoir densities which differ 
only slightly. Equation (f33|) leads to 



1 f Pl 
(j) = j dpD(p) 

L J PL 



(34) 



where p\ and pl are the densities at the first and 
the last site, respectively, and L is the length of the 
system. For \p\ — pj\ <C pi one has 



D 



Pl + PL 



L(j) 



(Pi - Pl)' 



(35) 
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FIG. 6: (a) Nonlinear diffusion coefficient D as a function of the density, obtained from the simulation data in the spreading 
regime (■) and in the steady state (□), as well as the corresponding analytical result [full line, Eq. ([56])]. All the data correspond 
to w = 1.5 and L = 1000. 

(b) Nonlinear diffusion coefficient [Eq. (|33[) ] from steady-state simulation data for w = 0.40 (+), w = 0.80 (□), and w = 1.00 
(■) (L = 1000). Analytical results (lines) for D(p) are calculated from Eq. (|56[) . 



B. Results from the spreading regime 

In order to extract D(p) from the numerical data for 
given w and L, we have considered all the profiles in 
the scaling regime. We have binned the p axis with 
Ap = 0.05, averaged all the values A in each bin, and 
evaluated the function X(p) by interpolation of the re- 
sulting data points, while jx/o(A) has been obtained by 
finite differences. The results for D(p) obtained by using 
Eq. are shown in Fig. EJa). While it appears that, 
for large values of p, D(p) — >■ 1 as expected from the 
corresponding discussion in Subsec. Ill Bl for p — ► the 
diffusion coefficient goes to zero due to the systematic er- 
ror in the derivative ^rp(A), as explained in Subsec. [IV Al 
The noise at large values of p is due to determining the 
derivatives numerically, because for large p the spatial 
fluctuations of the density are stronger. 

The diffusion coefficient is peaked and the substrate 
potential enhances diffusion (see Subsec. Hi Bp . The po- 
sition of the peak (p ~ 1) cannot be predicted by quali- 
tative arguments, but is in the range of low densities, as 
expected from the discussion in Subsec. Ill Bl 



tained from the spreading data and from the steady-state 
data is good for p > 0.5 and even better for p > 2, which 
shows that the diffusion picture described in Subsec. HV Al 
leads to consistent results. 

Simulations under quasi-equilibrium conditions have 
been restricted to the case w = 2 because the quantita- 
tive agreement between the results obtained from quasi- 
equilibrium, using Eq. (|35p for D(p), and those obtained 
in the steady state is satisfactory (see Fig. [7]). Due 
to the small difference in density of the two reservoirs 
(Sp = 0.01) the average current is very small and re- 
quires accurate measurements. The data are obtained 
by averaging over 10 7 configurations or more (approxi- 
mately 100 times more than for the steady state data), 
leading to a high precision for the density profile, too. 
The autocorrelation time for the average density has been 
carefully checked and the time intervals between samples 
have been chosen in order to minimize correlations. The 
procedure allows one to estimate reliably the statistical 
error for the diffusion coefficient, shown by the error- 
bars in Fig. [7fa). Note that the results obtained from 
steady-state simulations with bigger reservoir differences 
lie within the errorbars. 



C. Results from steady-state and quasi-equilibrium 
regimes 

In order to obtain D(p) from the steady-state data 
by using Eq. (|33j) we have measured the average cur- 
rent and the average density profile. The latter has 
been appropriately binned (the density is averaged over 
5 sites), in order to be able to evaluate the derivative via 
finite differences. The corresponding results are shown 
in Figs. Eta), [Spa), and[3 We note that D(p -> 0) -> 1. 
The correct behavior at low densities is captured, while 
for large p the data are still rather noisy, because the 
method to extract D relies on determining derivatives 
numerically. The overall agreement between results ob- 



V. CONTINUUM DESCRIPTION 

In this section we derive the nonlinear diffusion equa- 
tion corresponding to the continuum limit of the model. 
The equation is derived from the microscopic dynamics 
by using several simplifying assumptions. We start from 
the master equation, which describes exactly the dynam- 
ics of the model and in its most general form can be 
written as 

d t P t (C) = Y,M(C,C')P t (C), (36) 

c 

where C is a generic configuration, while A4(C, C) en- 
codes the transitions from the configuration C to the con- 
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FIG. 7: Nonlinear diffusion coefficient D(p) for large values of w. The open squares (□) are obtained from simulation data in 
the steady-state regime [Eq. (I33|l ] for w = 2 (a) and w — 4 (b) (L = 1000). The lines correspond to analytical calculations 
[Eq. (|56[l ]. The crosses (x) with errorbars in (a) are obtained from simulation data under quasi-equilibrium conditions [Eq. (|35l) ]. 



figuration C . In our case, the operator M can be split 
into bulk (Mb) and boundary (M s ) terms, so that 

d t p t (c) = \- M ^ c ') + Mb ^ c ')] p ^ c ')- ( 3? ) 

c 

The operator Mb describes bulk moves, upon which 
particles are exchanged between columns at sites x £ 
[1,...,L] and which are associated with the rates de- 
fined in Eq. (fT!?)) . The operator M s , describes boundary 
moves upon which particles are inserted into or removed 
from the system at the sites x — 1 or x — L, and which 
are associated with the rates introduced in Subscc. Ill CI 
Explicit expressions for the operators Mb and M s are 
given in Appendix |A"1 

The evolution of the ensemble average of a generic 
(time-indipcndent) operator O can be obtained from 
Eq. {31 as 



d t (0) = 0(C)M(C,C')P t (C), 



(38) 



where 0(C) is the value of the operator O for configura- 
tion C. Recalling that M(C,C) = - E C¥C M(C, C), 
it is straightforward to obtain 



d t {<D) = J2k(c)p(C), 

c 



(39) 



where IC(C) is the jump moment of the operator O de- 
fined as 

JC{C) = [0(C')-0(C)]M(C',C). (40) 

We consider now the operator Af x , defined as N X (C) — 
n x . Accordingly, the configurations C, for which the 
jump moments defined in Eq. (|4"0")) are nonzero, are C = 
{nx, • • • , n x -i + 1, n x — 1, n x +i, • • • , n L }, {rii, . . . , n x ^i - 



l,n x ,n x+ i + 1, ...,n L },{m,...,n x -i,n x + l,n x+ i - 
1, . . . , n L }, {nx, . . . , n x -i, n x - 1, n x+1 + 1, . . . , n L } with 
x £ [2, . . . , L - 1], so that using Eqs. (03), CHI), (|A"l~|) . 
and (IA2I) one obtains in the bulk 



d t p(x, t) = - [(j(x + 1, t) - j(x, t))] , (41) 

where p(x,t) = (n x ), x £ [2, ...,L — 1], and 

t)) = (^(na-i.na) - n x _i)) (42) 

is the mean local and instantaneous current in Eq. (|27[) . 
Assuming that the probability distribution [Eq. (pJT|) ] fac- 
torizes completely, i.e., within the mean field approxima- 
tion 



P(C,t) = Y[p y (n y ,t), 



(43) 



the average rates in Eq. ([42]) [for the definition of the 
rates see Eq. (fT!Zj) ] reduce to 



oo 



oo 



(u(n x ,n x+1 )) = X! "' 



Y\_p v (n y ,t) <j n K exp 

3/=l 



n i — 1 til-1 



(Tfe+i + l) 4 < 

fi(w,t,x)f 2 (w,t,x + 1), 



(44) 



where 



fx(w,t,x) = p a; (n 2; ,t)n a; exp 



and 



f 2 (w,t,x) = Y Px(n x ,t)exp 



n x =0 



w 



(n x + I) 4 



(45a) 



(45b) 
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The second equation in Eq. (|44|) holds because the sum 
in the first line equals 1 for any y ^ x, x + 1 due to 
the normalization of the distribution p. We further as- 
sume that the distribution p x depends smoothly on x 
and t, and that it does so only via the mean site den- 
sity p(x,t) (or, equivalently, on an effective local chem- 
ical potential) such that p x (n,t) — p(n,p(x,t)) and 
fi,2(w, t, x) — /i.2 {w, p(x, t)). Expanding the density up 
to second order in a [i.e., p(x+a) w p(x)+ad x p+\a 2 d 2 p] 
and then setting a = 1 [Eqs. (|TTj) and l|33])]. leads to the 
diffusion equation 

d t p(x,t) = [f 1 (w ) p)d p f 2 (w,p) - f 2 (w,p)d p f 1 (w,p)] 
x d*p(x,t), 

(46) 

with a density-dependent diffusion coefficient 

D(w,p) = fi(w,p)dpf 2 {w,p) - f 2 {w,p)d p fi(w,p). (47) 

In the steady-state regime the functions fi j2 can be com- 
puted explicitly in the local equilibrium approximation, 
i.e., by approximating the exact steady-state distribu- 
tion with a grand canonical equilibrium Gibbs distribu- 
tion Pq\ 



P G {n,w,p x 
where 



1 



1 



Z(w,p x ) n\ 



exp (2w h(n) — p x n) (48) 



Z(w, p x ) = V* — exp {2wh(n) - p x n); 



n=0 



(49) 



h is defined in Eq. (|16p and, as before, p x = (3p x is a 
dimensionless chemical potential. These approximations 
are expected to hold if the density varies slowly in space, 
so that the parts of the system to the left and to the 
right of x act on the column at x effectively as a particle 
reservoir with a well-defined chemical potential p x . Ac- 
cordingly, a nontrivial profile p x emerges which smoothly 
interpolates between the reservoir densities at x = and 
x = L + 1 

The distribution Pq can be expressed in terms of the 
local density p x by solving the implicit equation 



;M) = ^Pc{n,w,p x )n = p x 



(50) 



for p x = p(p x ). The functions f\ and f 2 in Eq. (I47|) are 
then given by 



1 



nl 



exp 



Z(w,p) 
2w h(n) — 



E 

n>0 



W 



(n + l) 4 



pn 



(51a) 



h{w,p) = r V 

Z(w,p)^ 

1 / , . W 

7^— TJj exp [2w h{n) + -j 
Due to Eq. fTB]) one has 



pn 



(51b) 



2w h(n) 



(n + l) 4 



which implies 



and thus 



2wh(n + 1) 



5u/ 2 = e-'U/j - / 2 . 



w 

(n + l) 4 ' 



(52) 



(53) 



(54) 



From Eq. (|50[) one can infer the derivative of p with re- 
spect to p: 



d^p = nd v. p G = E n W n ) ~ n ) Pg = 

n— 1 n — 1 

= (nf - (n 2 ) = -x{w,fx). 
Combining Eqs. (|4T|) . (f54"| . and ([55| one obtains 



(55) 



D(w, p) = 



fi(w,p)f 2 (w,p) 



X(w,p) 
(u (n s ,n x+ i)) 



m=m(p) 



X{w,p) 



(56) 



Equation (156p allows one to compute numerically the dif- 
fusion coefficient. To this end we solve Eq. ([50)1 for p(p) 
and insert the solution p(p(x)) back into Eq. ([56]) . The 
resulting Z?(p) is obtained by approximating the series 
in Eq. (|5"T|) by finite sums. We have checked the stabil- 
ity of the calculation for densities < p < 11, in order 
to be able to make contact with our Monte Carlo data. 
For w < 1.5 the theoretical expression is in good agree- 
ment with the simulation data (see Fig. [5]), but for larger 
values of the interaction significant deviations occur (see 
Fig. [7|). These deviations systematically increase with 
increasing interaction strength, which cannot be easily 
blamed on numerical inaccuracies. A first possible ex- 
planation for the deviations could be the non-equilibrium 
character of the simulations, due to the chemical poten- 
tial gradient present in the system. However, the quasi- 
equilibrium simulation results allow us to rule out strong 
non-equilibrium effects as the primary source for these 
deviations, because D(p) computed from these data co- 
incides with the one obtained from non-equilibrium sim- 
ulations; thus, the diffusion coefficient is basically inde- 
pendent of the difference between the reservoirs densi- 
ties. Taking advantage of this fact we can make use 
of general results for an infinitely large system, derived 
in quasi-equilibrium conditions, such as the Green-Kubo 



12 



formula 142 



D(p) 



{u(n x , n x+ l))eq 



°o poo 

Y] / dt' (j x+1 U eq (t')j x :) eq 
x'=l Jo 



(57) 



where {-) eq indicates the average performed over the equi- 
librium distribution and U eq is the evolution operator for 
the equilibrium dynamics. Following the discussion in 
Ref. [42|, in Appendix lAl we present a brief derivation of 
Eq. (|57[) for the class of models we are interested in. The 
comparison of Eqs. (|55|) and ([57)) shows that our mean- 
field calculation reproduces the first term in Eq. (|57|) . 
while the terms which would reduce (in the appropriate 
limit) to the time integral of the current correlations can- 
not be captured by this mean-field approximation. 

Computing explicitly the current correlations is dif- 
ficult, but some of their general properties [HJ allow 
us to conclude that they provide a qualitatively cor- 
rect correction to the diffusion coefficient calculated via 
Eq. (151)1) . The function (j x +iU e q(t)j x >) e q appearing in 
Eq. (|57p is integrable, positive, and decaying exponen- 
tially for t — > oo. Thus it is a negative contribution to 
Eq. (|57[) and decreases the diffusion coefficient obtained 
from Eq. ([55)) . This is in qualitative agreement with the 
data shown in Fig. [7] 



VI. 



SUMMARY AND CONCLUSIONS 



We have introduced a lattice model (Fig. [J) for spread- 
ing of a fluid in narrow, quasi one-dimensional slit-like 
channels and in contact with particle reservoirs located 
at their ends. The model accounts for long-range attrac- 
tive substrate-fluid interactions, while the fluid-fluid in- 
teraction is taken to be hard-core only. We have studied 
the spreading behavior and stationary state using kinetic 
Monte Carlo simulations and a non-linear diffusion equa- 
tion corresponding to the continuum limit of our discrete 
model. The main results are the following. 

The spreading regime has been studied starting from 
an empty lattice: we have set the right reservoir density 
to pl+i = 0, so that it acts as a perfect sink, and the left 
reservoir to a nonzero value (typically po — H)j thereby 
feeding particles into the system. At intermediate times, 
for which the reservoirs do not play a relevant role, we 
have found a diffusion-like behavior, such that profiles at 
different times collapse onto a single master curve if the 
scaling variable A = xj \ft is introduced (see Fig. [3|) . This 
scaling is further confirmed by the time evolution of the 
total mass [see Eq. ([25]) ] shown in Fig. [4] In the steady 
state (Fig. [5]) we have found nontrivial density profiles 
depending on the substrate interaction strength. The 
dynamics at the boundary influences the density profiles: 



the analysis in Subsec. lll Cl shows that one of the possible 
definitions of the boundary rates causes discontinuities in 
the profiles. These discontinuities can be eliminated by 
an alternative definition or can be calculated explicitly 
(see Fig. [2]). The profiles also depend on the reservoir 
density as discussed in Appendix [Bl either lying com- 
pletely above the free-diffusion straight line or crossing 
this line at a certain point. This feature can be described 
by the deviation Ax(p) from free diffusion [see Eq. (|B1[) 
and Fig. [5] and explained in terms of the function R(p) 
introduced in Eq. (|B8[) . In Fig. [TU] the general behavior 
of this function is sketched while in Fig. [5] the simulation 
and mean-field results for the deviation Ax(p) are shown. 

Assuming that a nonlinear diffusion equation describes 
the behavior of the particle density, we have extracted 
the diffusion coefficient from the simulation data in the 
spreading regime (Fig. [UJa)) and in the steady state 
(Figs. [H] and [7]) . The two sets of results are compatible 
with each other and show that the interaction with the 
substrate tends to enhance diffusion, as expected from 
qualitative arguments discussed in Subsec. Ill Bl 

The Monte Carlo results for the diffusion coefficient are 
in agreement with an analytical calculation based on local 
equilibrium assumptions, in which the equilibrium grand 
canonical distribution is modified by a spatially varying 
local chemical potential. The agreement is very good for 
weak interactions with the substrate (Fig. [5]), but deteri- 
orates for strong interactions (Fig.[7J). The independence 
of the diffusion coefficient from the difference between 
the reservoirs densities, suggested by quasi-equilibrium 
simulations, allows one to qualitatively explain these de- 
viations by using a Green-Kubo formula. 

An experimental situation to which our model would 
apply should be effectively quasi one-dimensional, with 
effective inert side walls confining the fluid and the 
substrate-fluid interaction strongly dominating over the 
fluid- fluid interaction. This might be realized, e.g., by 
colloidal particles in channels. An experimental investi- 
gation of the diffusion coefficient as a function of the den- 
sity would allow one to make a direct comparison with 
the results presented here. 



APPENDIX A: A GREEN-KUBO FORMULA 



In this appendix we derive a Green-Kubo formula in 
the case of the non-equilibrium steady state induced by 
an infinitesimal difference 5^ = po — pl+i — ► between 
the dimensionless (in units of fc^T) chemical potentials 
po and pl+i of the left and the right reservoir, respec- 
tively. The rates are denoted by u(n x , n x +i) but they are 
not necessarily of the explicit form of Eq. (fT2|) . The mas- 
ter equation, which describes the dynamics of the system, 
is given by Eq. (|37p. where M. & is an operator that acts 
on a distribution P as 
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L-l 

5^ M b (C, C")P(C") = J2 {- I w n -+i) + M K+i. n -)l p ("i. • • • . u l) ■ 

C x=l 

+u (n x+ i + l,n x - l)P(m, . . .,7ia: - 1,^+1 + 1, ■■■,n L ) 
+u (n x + l,n x+ i - l)P(m, . . .,n x + l,n x+ i - I,.. .,n L )} , 



r 



(Al) 



and where M s acts on P as 



M S (C,C')P(C) 



- [u (ni,n ) +u(n ,ni) + u (ni,n i+ i) + u (ul+u n L )] P (m, ■■■,n L ) 
+u (no, ni - 1) P (ni - 1, . . . , n L ) + u (nj + 1, n ) P (ni + 1, . . . , n L ) 
+it (n_L + 1, n i+ i) P (m, . . . , m, + 1) + u (n L +i,n L - 1) P (ni, . . . , n L - 1). 



r 



(A2) 



If the alternative definitions of the rates at the bound- 
aries are employed [Eq. pijl] the rates u ayK . n j replace 
u(n ,ni) ,u(ni,n ) , and u(n Ll n L+1 ) ,u(n L+1 ,n L ) in 
Eq. (32). 

We require that the rates u satisfy detailed balance, 
so that for 5^ = 0, i.e., in the equilibrium case and thus 
without net transport of particles, the system is described 
by the grand canonical distribution 



pGC 



pC -fi N 
1 eq c > 



where P^ q is a canonical distribution and N = ^2 X=1 n x 
is the total number of particles in the system. In the 
general case 8^ ^ detailed balance should be satisfied 
locally by the boundary rates at each end of the system. 
The formal solution of the master equation (|37|) is 



P(t) 



P. 

J- II 



U{t)P in , 



(A3) 



where Pi n is the initial distribution P(t = 0) and U(t) = 
e Mt is the time evolution operator. In the equilibrium 
case 8^ = the corresponding dynamics is indicated by 
M eq . Splitting arbitrarily the operator M as M — Mq + 
Mi, one can check that from the equation 

U(t) = U (t) + f dt'U(t - t')MiU Q {t'), (A4) 
Jo 

where Uo(t) = e Mot , the two equations 
d t U(t) = MU(t) 
17(0) = 1 



(A5) 



can be obtained. Since these equations have as unique 



solution the evolution operator U(t) = e M , Eq. (|A4[) 
holds. Assuming that unique stationary distributions Pq 
and P st exist for the dynamics defined by Mo and M, 
respectively, i.e., 



M P = 0, 
MP st = 0, 



(A6) 



Eqs. (3I|) and (36j) lead to 

P st = lim [U(t)P ] = P + / dtU(t)MP . (A7) 

The arbitrariness in the splitting of M allows one to 
choose Pq as the local equilibrium distribution 



P a = P^ q exp 



x=l 



(A8) 



For ^ < 1 one has fj,(x) ~ \i + S^x/L); applying M 
to P [sec Eq. (jAll) ] and by using the detailed balance 
condition one obtains 



L-l 



MP = ^ [ u ( n x,n xA 



e L — 1 
e T — 1 



+u(n x+1 ,n x ) ei -1 j Po({ni,n L }) 



(A9) 

The aim is to calculate the mean current (j) s t, from 
which one can obtain D in the limit of small density 
differences by using Eq. (|35[) . Averaging j, given by 
Eq. (|42f , over the distribution P st and expanding to first 
order in 8 leads to 



(J 



x+l/st 



(u(n x ,n x+ i))o 




(A10) 



where (-) indicate the average taken over the local equi- 
librium distribution [Eq. (|A8[) ]. By taking the limits 
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Sfj, — > and L 



oo we obtain the Green-Kubo formula 



D(p) 



(u(n x ,n x+ i)) 



eq 



oo poo 

/ dt' (j x+1 U eq (t')j x ,) t 

x'=l J ° 



(All) 



The limits <5 M — > and i 



withx = (n x ) e q-{n x )l q . 

imply t hat P -> P eg and U -> L/ eg = e^ 9 *. Note that 
m Eq. (|ATT|) the averages are taken over the equilibrium 
distribution and the properties of the diffusion coefficient 
are completely determined by the equilibrium dynamics 
of the model. Accordingly, the r.h.s. of Eq. (|A11[) is 
independent of x and thus D depends on p only. 



APPENDIX B: INFLUENCE OF THE 
RESERVOIR DENSITIES ON THE DENSITY 
PROFILES 




FIG. 8: Schematic comparison between a typical steady-state 
density profile (curved line) and free diffusion (straight line). 
At a given density p the corresponding positions x(p) and 
xo(p) are taken from the steady-state density profile and from 
the free-diffusion line, respectively, and Ax(p) is determined 
from Eq. ilBil . 



In Sec. IIVI we have shown that, within our model, the 
density-dependent diffusion coefficient D(p) does not de- 
pend on the boundary conditions, i.e., the densities of the 
left and right reservoirs. As shown in Fig. [5l the density 
profiles for low reservoir densities seemingly lie mostly 
above the free diffusion straight line [see Fig. [5jb)] , 
whereas the opposite happens for high reservoir densi- 
ties [see Fig. [SJa)]. While these profiles refer to different 
interaction strength w, we note that the difference be- 
tween D(p) for w = 2 and for w = 1.5 [see Figs. HJb) 



and [TJa)] is not very marked. Therefore it is unlikely 
that this difference explains the differences between the 
profile corresponding to w = 1.5 in Fig. [5ja) and the 
profile corresponding to w = 2 in Fig.[5]Jb). 

In order to clarify the relation between the steady-state 
profile and the boundary conditions, here we restrict our 
analysis to the case in which the right reservoir is a per- 
fect sink, (pl+i — 0) and the left reservoir density po is 
varied, keeping the interaction strength w constant. The 
definitions in Eq. ([21) for the rates at the boundaries 
have been used, so that the density profile is a continu- 
ous function of x. 



1 



X 

< 

X 

o 



0.5 







-0.5 




FIG. 9: Ax(p) [see Eq. (|B1|) ] for left reservoir densities po — 1 
(+), p = 2.5 (□) and p = 4 (0) with fixed w = 1.5, L = 300, 
and pL+i = 0. The full lines correspond to Ax(p) obtained 
by using the mean-field diffusion coefficient [Eq. (|56| . The 
quantity Ax(p) in the y-axis is multiplied by a factor 10 in 
order to enhance visibility. 



A convenient quantity to describe the dependence on 
po of the profiles is the deviation in x at fixed p be- 
tween the actual density profile and the straight line cor- 
responding to free diffusion (see Fig. [5]), defined as 



Ax(p) 



x(p) - x (p) 
L + l ' 



(Bl) 



where x(p) is the inverted steady-state profile and 
xo(p) = (L + 1)(1 — p/ po) is the corresponding free 
diffusion line. In Fig. [S] we report the numerical data 
for Ax{p) corresponding to three left reservoir densities 
(po = 1, 2.5, and 4). These results exhibit the typical be- 
havior of the deviation Ax: at small po it is positive for all 
p (po = 1) whereas upon increasing the reservoir density 
a negative part appears (po — 2.5 and 4). This pattern 
can be explained along the line of arguments introduced 
in Sec. IIVI based on the coarse-grained description of the 
diffusion-like behavior. The nonlinear diffusion equation 
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[Eq. ([30f ] yields for the inverse steady-state profile x(p) 
dx(p) _ D(P) 



dp 



J 



(B2) 



where j is the current and D the diffusion coefficient. 
Integrating Eq. (|B2|1 one obtains 



1 f Po 
x{p) = - / d£D(0, 
J Jp 



(B3) 



where < p < po- Imposing the boundary condition 
x(Q) = L + 1 at the right end of the system leads to 



1 



L + l 



(B4) 



so that x(p) is given by 



x(p) = (L + l) 



d^(0 



d£D(0. (B5) 



By combining Eqs. (|B1[) and (|B5P one obtains 
z(p) - ^o(p) 



Aar(p) 



L+l 



Po - P 
Po 



(B6) 



which, as expected, vanishes for p = and p = pq. One 
can re- write the last equation as follows: 



Ax(p) = p 



I'D 



i r Po l r p 

- / d£D £ - - / d£D(£) 
Po 7 P Jo 



(B7) 



Note that p [J Q pa .D(£)d£] 1 is nonnegative whereas the 
second factor on the r.h.s. of Eq. (|B7|) can change sign. 
Denoting this latter factor by h{p) and introducing the 
function f(p) — D(p) — 1, one obtains 

Hp) = - r d^/(0 - - = R(po) - R(p), 

Po Jo P Jo 

(B8) 



where R{p) = ± J Q P d£ - 1]. Within our model the 

function / is always positive (apart from Fig. [SJa) for 
p — > and some noisy data at large p in Figs. [6] and [7]) 
and it vanishes for p — > and p — > oo. Thus also 
vanishes for p — > or p — > oo and hence has a maximum 
at a certain pm, with < pm < oo. A sketch of this 
function is provided in Fig. 1101 Two qualitatively differ- 
ent behaviors emerge by varying the left reservoir density 
Pq. If the reservoir density is set to p^ < p^, then in the 
interval < p < p^ no solution to R(pq) — R(p) — can 
be found other than the trivial one p = Pq. In this situ- 
ation, we have R(p) < R{Pq ) so that the function Ax(p) 
is always positive in the considered interval of densities, 




R{p 



FIG. 10: Qualitative plot of the function R(p) (see Eq. jBgl ). 
Two possible situations are shown: if the reservoir density 
Pq is smaller than pM, there are no densities < p < p$ 
such that R(p) = R(pq). In the opposite case with reservoir 
density pg" > pm, there exist a density < ps < Pq such that 
R( PS )=R(p>). 



and the density profile always lies above the line corre- 
sponding to free diffusion. This is the case for po — 1 and 
substrate interaction w = 1.5 as shown in Fig. [5J In the 
opposite situation, in which the density of the reservoir 
is set to Pq > pm , there always exists a density ps < Pq 
such that R(pq ) — R{ps) = 0. At this density the den- 
sity profile crosses the free diffusion line. For p > pg, 
Ax(p) > while Ax(p) < for p < p s . In Fig. [5] this 
corresponds to the cases po = 2.5 and po = 4. 
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